%% 参数设置
fs = 10000;       % 采样率 (Hz) - 需满足奈奎斯特（>2 * 10kHz）
t = 0:1/fs:1-1/fs; % 时间向量 (1秒)
f_excite = 10000; % 激励频率 (Hz)
f_interf = 50;    % 干扰频率 (Hz) - 如工频
A = 1;            % 激励幅值
B = 0.2;          % 干扰幅值
phi = pi/4;       % 干扰相位偏移

%% 生成信号
t_excite = 0:1/fs:1-1/fs; % 时间向量
s = A*sin(2*pi*f_excite*t_excite); % 激励信号
n = B*sin(2*pi*f_interf*t_excite + phi); % 干扰信号
x = s + n; % 含干扰的测量信号（可添加噪声：x = x + 0.05*randn(size(t_excite))）

%% 傅里叶变换分析
N = length(x);
X = fft(x);
f = (0:N-1)*(fs/N); % 频率轴

%% 干扰频谱识别
[~, idx] = findpeaks(abs(X), 'SortStr', 'descend', 'NPeaks', 1);
interf_freq = f(idx); % 提取干扰频率（可能需调整阈值）

%% 生成反向干扰信号
phase_interf = angle(X(idx)); % 提取干扰相位
n_reversed = -B*sin(2*pi*interf_freq*t_excite + phase_interf); % 精确反向

%% 干扰抵消
x_clean = x + n_reversed;


%% 频域验证
X_clean = fft(x_clean);
figure;
subplot(3,1,1);
plot(t_excite, x); title('原始含干扰信号');
xlabel('时间 (s)'); ylabel('幅值');

subplot(3,1,2);
plot(f, abs(X)); title('原始信号频谱');
hold on;
plot(interf_freq, abs(X(idx)), 'ro', 'MarkerFaceColor', 'r');
xlim([0 60]); % 聚焦低频干扰
xlabel('频率 (Hz)'); ylabel('幅值');

subplot(3,1,3);
plot(t_excite, x_clean); title('抵消后信号');
xlabel('时间 (s)'); ylabel('幅值');